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Total Internal Reflection Fluorescence Cross Correlation Spectroscopy (TIR-FCCS) has recently 
(S. Yordanov et al., Optics Express 17, 21149 (2009)) been established as an experimental method to 
probe hydrodynamic flows near surfaces, on length scales of tens of nanometers. Its main advantage 
is that fluorescence only occurs for tracer particles close to the surface, thus resulting in high sensitiv- 
ity. However, the measured correlation functions only provide rather indirect information about the 
flow parameters of interest, such as the shear rate and the slip length. In the present paper, we show 
how to combine detailed and fairly realistic theoretical modeling of the phenomena by Brownian 
Dynamics simulations with accurate measurements of the correlation functions, in order to establish 
a quantitative method to retrieve the flow properties from the experiments. Firstly, Brownian Dy- 
namics is used to sample highly accurate correlation functions for a fixed set of model parameters. 
Secondly, these parameters are varied systematically by means of an importance-sampling Monte 
Carlo procedure in order to fit the experiments. This provides the optimum parameter values to- 
gether with their statistical error bars. The approach is well suited for massively parallel computers, 
which allows us to do the data analysis within moderate computing times. The method is applied 
to flow near a hydrophilic surface, where the slip length is observed to be smaller than lOnm, and, 
within the limitations of the experiments and the model, indistinguishable from zero. 
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I. INTRODUCTION 



A good understanding of liquid flow in confined ge- 
ometries is not only of fundamental interest, but also im- 
portant for a number of industrial and technological pro- 
cesses, such as flow in porous media, electro-osmotic flow, 
particle aggregation or sedimentation, extrusion and lu- 
brication. It is also essential for the design of micro- and 
nano-fluidic devices, e. g. in lab-on-a-chip applications. 
However, in all these cases, an accurate quantitative de- 
scription is only possible if the flow at the interface be- 
tween the liquid and the solid is thoroughly understood 
p]-[l0| . While for many years the so-called no-slip bound- 
ary condition (relative velocity at the interface equal to 
zero) had been successfully applied to describe macro- 
scopic flows, more recent investigations revealed that this 
condition is insufficient to describe the physics when fiows 
through channels with micro- and nano-sizes are consid- 
ered [3, Q. On such small scales, the relative contri- 
bution from a residual slip between liquid and solid be- 
comes important. This is commonly described by the 
so-called slip boundary condition, which is characterized 
by a non- vanishing slip length Is, defined as the ratio of 
the liquid dynamic viscosity and the friction coefficient 
per unit area at the surface. An equivalent definition is 
obtained by taking the ratio of the finite surface fiow ve- 
locity, the so-called slip velocity Vg , and the shear rate at 
the surface: Is = Vs/{dv/dz)z=o, where z is the spatial 
direction perpendicular to the surface, located at z = 0. 
This boundary condition is the most general one that is 



possible within the framework of standard hydrodynam- 
ics [ll|; the no-slip condition is simply the special case 

ls=0. 

The experimental determination of the slip length, 
however, is challenging, since high resolution techniques 
are needed to gain sufficiently accurate information close 
to the interface. Hence, the existence and the magnitude 
of slip in real physical systems, as well as its possible de- 
pendence on the surface properties, are highly debated 
in the community, and no consensus has been reached so 
far. Clearly, a resolution of these controversies requires 
further improvement of the experimental techniques. 

To date, two major types of experimental methods, of- 
ten called direct and indirect, have been applied to study 
boundary slip phenomena. In the indirect approach, an 
atomic force microscope or a surface force apparatus is 
used to record the hydrodynamic drainage force neces- 
sary to push a micron-sized colloidal particle versus a 
flat surface as a function of their separation [l^l ■ The 
separation can be measured with sub-nanometric resolu- 
tion, and the force with a resolution in the pN range. A 
high force is necessary to squeeze the liquid out of the 
gap if the mobility of the liquid is small. Conversely, if 
the liquid close to the surface can easily slip on it, then a 
small force is necessary. From this empirical observation 
a quantitative value of the slip length can be deduced 
using an appropriate theoretical model [1, 0, . While 
this approach is extremely accurate at the nanoscale, it 
does not measure the flow profile directly. 

Direct experimental approaches to flow proflling in mi- 
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crochannels are commonly based on various optical meth- 
ods to monitor fluorescent tracers moving with the liquid. 
Basically they can be divided into two sub-categories. 

The imaging-based methods use high-resolution opti- 
cal microscopes and sensitive cameras to track the move- 
ment of individual tracer particles via a series of images 
P"5l - l2ll | . While providing a real "picture" of the flow, the 
imaging methods have also some disadvantages related 
mainly to the limited speed and sensitivity of the cam- 
eras: relatively big tracers are needed, the statistics is 
rather poor, and large tracer velocities cannot be easily 
measured. 

In the Fluorescence Correlation Spectroscopy (FCS) 
based methods the fluctuations of the fluorescent light 
emitted by tracers passing through a small observation 
volume (typically the focus of a confocal microscope) are 
measured [22|. Using correlation analysis and an appro- 
priate mathematical model the tracers' diffusion coeffi- 
cient and flow velocity can be evaluated [23l - [26j . In par- 
ticular, the so-called Double-Focus Fluorescence Cross- 
Correlation Spectroscopy (DF-FCCS) that employs two 
observation volumes (laterally shifted in flow direction) 
is a powerful tool for flow profiling in microchannels p?! - 
Isoj . Due to the high sensitivity and speed of the used 
photo detectors (typically avalanche photodiodes) in the 
FCS based methods even single molecules can be used 
as tracers. Furthermore, the evaluation of the velocity is 
based on large statistics and high tracer velocities can be 
measured. 

During the last two decades both the imaging and the 
FCS methods were well developed to the current state 
that allows fast and accurate measurements of flow ve- 
locity profiles in microchannels. The situation, however, 
is different when the issue of boundary slip is consid- 
ered. Due to the limited optical resolution imposed by 
the diffraction limit, it is commonly believed that these 
methods are less accurate than the force methods dis- 
cussed above and cannot detect a slip length in the tens of 
nanometers range. On the other hand, the possibility to 
directly visualize the flow makes the optical methods still 
attractive and thus continuous efforts were undertaken 
to improve their resolution. One of the most successful 
approaches in this endeavor is Total Internal Reflection 
Microscopy (TIRM) [3l| , which uses total internal reflec- 
tion at the interface between two media with different 
refractive indices, like, e. g., glass and water. This cre- 
ates an evanescent wave that extends into the liquid in 
a tunable region of less than ~ 200nm from the inter- 
face. Optical excitation of the fluorescent tracers is then 
possible only within this narrow region. During the last 
few years TIRM was successfully applied to improve the 
resolution of particle imaging velocimetry close to liquid- 
solid interfaces [l8l - l2]| . and slip lengths in the order of 
tens of nanometers were evaluated. With respect to FCS, 
however, TIR illumination had, until recently, been lim- 
ited to diffusion studies only [32|, [s^ , while there were no 
reports on TIR-FCS based velocimetry and slip length 
measurements. 



With this in mind, we have recently developed a new 
experimental setup that combines for the first time TIR 
illumination with DF-FCCS for monitoring a liquid fiow 
in the close proximity of a solid surface [s^. Such a 
combination offers high normal resolution, extreme sen- 
sitivity (down to single molecules), good statistics within 
relatively short measurement times and the possibility 
to study fast fiows. Our preliminary studies have shown, 
however, that the accurate quantitative evaluation of the 
experimental data obtained with this TIR-FCCS setup is 
not straightforward because the model functions needed 
to fit the auto- and cross-correlation curves (and extract 
the fiow velocity profile) are not readily available. The 
standard analytical procedure to derive these functions is 
[27l - l29l |: (i) solve the convection-diffusion equation with 
respect to the concentration correlation function, (ii) in- 
sert the derived solution in the corresponding correlation 
integral and (iii) solve it to finally get the explicit form 
of the correlation functions. This pro cedure was success- 
fully used by Brinkmeier et al. [231 to derive analyti- 
cal expressions for the auto- and cross-correlation func- 
tions obtained with double focus confocal FCCS (i. e. 
with focused laser beam illumination as opposed to the 
evanescent wave illumination in our case), where it was 
assumed that the flow velocity and tracer concentration 
are spatially constant, which simplifles the calculation 
substantially. Such an assumption is reasonable if the 
observation volumes (laser foci) are far away from the 
channel walls, in the same distance. In the case of TIR- 
FCCS, however, the situation is different: The experi- 
ments are performed in the proximity of the channel wall 
and the distribution of the fiow velocity inside the ob- 
servation volume has to be considered. Furthermore, the 
concentration of tracers may also depend on z due to 
electrostatic repulsion or hydrodynamic effects. Finally 
the presence of a boundary, which must also be taken 
into account in the theoretical treatment, further com- 
plicates the problem. Therefore, a faithful description of 
the physics of TIR-FCCS makes the problem of calculat- 
ing the correlation functions (rather likely) unsolvable in 
terms of closed analytical expressions. 

For this reason, we rather resort to numerical methods, 
and in the present paper describe and test the procedure 
that we have developed: We employ Brownian Dynam- 
ics techniques to simulate the tracers' motion through 
the observation volumes and generate "numerical" auto- 
and cross-correlation curves that are consequently used 
to fit the corresponding experimental data. This fitting is 
done via Monte Carlo importance sampling in parameter 
space. The method is therefore fully quantitative, while 
not being hampered by any difficulties in doing analytical 
calculations. It should be noted that this approach also 
provides a substantial amount of flexibility: The details 
of the physical model are all encoded in the Brownian 
Dynamics simulation which specifies how the tracer par- 
ticles move within the flow. In the present work we have 
assumed a simple Couette fiow with a finite slip length, 
while the particles are described as simple hard spheres 
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with no rotational degree of freedom, and no interaction 
with the wall except impenetrability. It is fairly straight- 
forward to improve on these limitations, by, e. g., includ- 
ing hydrodynamic and electrostatic interactions with the 
wall, rotational motion of the spheres, or polydispersity 
in the particle size distribution. Moreover, the geom- 
etry of the observation volumes can be easily changed 
as well, and we have made use of this possibility in our 
present work, but only to some extent. Further refine- 
ments are left for future work, in which the basic method- 
ology would however remain unchanged. 

To test the accuracy of the newly developed TIR-FCCS 
experimental setup and the numerical data evaluation 
procedure, we have studied aqueous flow near a smooth 
hydrophilic surface and evaluated the slip length to be 
between and lOnm (however with a systematic error 
that is hard to quantify, and whose elimination would 
need a more sophisticated theoretical model). It is com- 
monly accepted \XL,1S. 21, 35 38| that the boundary slip 
should be zero (or very small) in this situation. Thus, our 
results indicate that TIR-FCCS offers unprecedented ac- 
curacy in the lOnm range for the measurement of slip 
lengths by an optical method. We believe that our re- 
sult for the slip length will be fairly robust, even if the 
physical model is refined further. 

Section |ll] outlines the experimental setup, while Sec. 
mil presents the experimental results and the numerical 
fits. We find that the measured cross-correlation func- 
tions deviate considerably from the model functions at 
short times, probably as a result of some optical effects 
which at present we do not fully understand. However, 
we show a practical way to eliminate such effects to a 
large extent, by means of a simple subtraction scheme. 
The following parts then outline in detail how the theo- 
retical curves have been obtained: Firstly, Sec. |IV] eluci- 
dates the relation between the measured correlation func- 
tions and the underlying dynamics of the tracer particles. 
We then proceed to describe the Brownian Dynamics al- 
gorithm to sample the model correlation functions (Sec. 
W} . Section lVTl then provides a detailed theoretical analy- 
sis of our subtraction scheme. In Sec 



Vlll we describe the 
Monte Carlo method to find optimized parameter values 
of our model. Section [Villi then discusses our results, in 
particular concerning the slip length; this is followed by 
a brief summary of our conclusions (Sec. IIX|) . 



II. EXPERIMENTAL SETUP 

Since the TIR-FCCS experimental setup has already 
been described in great detail elsewhere |34l | , only a brief 
qualitative overview of the basic ideas and quantities 
is given below. A scheme of the experimental setup is 
shown in Fig. [T] It is based on a commercial device 
(Carl Zeiss, Jena, Germany) that consists of the FCS 
module ConfoCor2 and an inverted microscope Axiovert 
200. The TIR illumination is achieved by focusing the ex- 
citation laser beam (488nm, Ar"*" Laser) on the periphery 




FIG. 1: (Color) Scheme of the experimental TIR-FCCS setup. 
BFP - back focal plane of the objective; DM - dichroic mir- 
ror; M50/50 - neutral 50% beam splitter; EFl, EF2 - emis- 
sion filters; PHI, PH2 - pinholes; APDl, APD2 - avalanche 
photodiodes; LI - tube lens; L2 - collimator lens; M - collima- 
tor's prism based mirror. Note that the two spatially sepa- 
rated observation volumes are created by shifting the pinholes 
PH1/PH2 in the x-y-plane. The cyan color indicates the exci- 
tation wavelength and the yellow-green color the fluorescence 
light, respectively. 



of the back focal plane (BFP) of an oil immersion mi- 
croscope objective with numerical aperture NA = 1.46. 
This leads to a parallel laser beam which emerges out 
of the objective and then enters the rectangular flow 
channel through its bottom wall (Fig. [T]). By adjust- 
ing the angle of incidence above the critical angle (« 61° 
for the glass-water interface) total internal reflection is 
achieved. In this situation only an evanescent wave ex- 
tends into the liquid and can excite the fluorescent tracers 
suspended in it. The intensity distribution of this wave 
in the x-y-plane (parallel to the interface) is Gaussian 
with a diameter of ^ 30/iTO (at e~^). In the z direction 
the intensity decays exponentially, I{z) — Iq cxjp{— z / dp) . 
The characteristic decay length dp, also called penetra- 
tion depth, depends on the laser wavelength A, the re- 
fraction indices of both media (rii - glass, n2 - water) 
and can be varied in the range 80 — 200nm by changing 
the angle of incidence. Thus the evanescent wave can 
excite only the tracers flowing in the proximity of the 
channel wall. The produced fluorescence light is collected 
by the same microscope objective and is equally divided 
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FIG. 2: (Color) The coordinate system and the linear flow 
field employed in the TIR-FCCS experiment. Wi and W2 
denote the shape and location of the observation volumes as 
seen by pinhole PHI and pinhole PH2, respectively; dp is the 
penetration depth which defines the axial extent of the obser- 
vation volume; wq is the typical extension of the observation 
volumes in the a;-y-plane; s^; indicates the observation vol- 
umes separation, center-to-center distance; v^c is the velocity 
field in positive x direction, which depends linearly on z. 



by passing through a neutral 50% beam splitter to enter 
two independent detection channels. In each channel the 
fluorescent light passes through an emission filter and a 
confocal pinhole to finally reach the detectors, two single 
photon counting avalanche photodiodcs (APDl, APD2). 
The pinholes PHI and PH2 define two observation vol- 
umes that are laterally shifted with respect to each other 
along the flow direction as schematically shown in Fig. 
[5J The center-to-center distance Sx between the two ob- 
servation volumes can be continuously tuned from to 
The signals from both channels are recorded and 
correlated to flnally yield the auto- and cross-correlation 
curves that contain the entire information about the flow 
properties, slip length and shear rate, close to the inter- 
face. 

The experiments were performed with a rectangu- 
lar microchannel of Ly = imm width, ~ lOOfim 
height and = 50mm length fabricated using a three- 
layer sandwich construction as described in earlier work 
[11 nil. The bottom channel wall at which the TIR- 
FCCS experiments were performed was a microscope 
cover slide made of borosilicate glass with a thickness of 
llOfim, cleaned with 2% aqueous solution of HcUmanex 
and Argon plasma. The root-mean-square roughness of 
the glass surface was in the range of 0.3nm and the water 
advancing contact angle below 5° (hydrophilic surface). 
The flow was induced by a hydrostatic pressure gradi- 
ent, created by two beakers of different heights, where 
the water level difference was kept constant by a pump. 
This allowed us to vary the shear rate near the wall in 
the range — 5000s~"'^. 



Carboxylate-modiflcd quantum dots (Qdot585 ITK 
Carboxyl, Molecular Probes, Inc.), with a hydrodynamic 
radius Rh = 6.87nm, were used as fluorescent trac- 
ers. The particles were suspended in an aqueous solution 
of potassium phosphate {K2HPO4) buffer [pH ~ 8.0, 
concentration 6mA/). The concentration of the quan- 
tum dots was found from our data analysis (see below) 
as ^ 30nM, corresponding to roughly 18 particles per 
(/irn)"^. 



III. CORRELATION CURVES 

The motion of the fluorescence tracers results in two 
time-resolved fluorescence intensities /i(t) and l2{t), 
which were measured with the two photo detectors. For 
the present system, we may safely assume that it is cr- 
godic and strictly stationary on the time scale of the 
experiment, such that only time differences matter [39j . 
Therefore, we may deflne the intensity fluctuations via 

5h{t) = im~{h), (1) 

where (•) denotes a time average or, equivalently, an en- 
semble average, and evaluate the time-dependent auto- 
and cross-correlation functions via the definition 

It should be noted that possible small differences in the 
sensitivity of the photo detectors or in the illumination 
of the pinholes cancel out, since in Eq. [5] only ratios of 
intensities occur. Gn and G22 are the two autocorrela- 
tion functions of pinholes 1 and 2, respectively, while 6*12 
and G21 are the forward and backward cross-correlation 
functions, respectively. It should be noted that in the 
presence of flow G12 and G21 differ substantially. In the 
limit of the two pinholes being located at the same posi- 
tion, the intensities Ii and I2 coincide, such that in this 
case all four entries of the matrix Gij are identical. 

Figure [3] summarizes our experimental results for the 
Gij and / or linear combinations thereof. Concerning 
the autocorrelation curves Gn and G22, we find that 
they are practically identical, which means that for the 
modeling it is safe to assume that both pinholes have 
the same properties. This is clearly shown in part (a), 
where one sees that Gn — G22 differs only marginally 
from zero (while in our model we have anyway strictly 
Gil = G22). Therefore, wc just used the arithmetic mean 
(Gil + G22)/2 (part (b)) as autocorrelation input for our 
fits, while wc discarded the Gn — G22 data. Concerning 
the cross-correlations, one sees that the forward function 
G12 (see part (c)) exhibits a pronounced peak, which is 
indicative of the typical time that a particle needs to 
travel from observation volume 1 to observation volume 
2. Another striking feature of G12 is the large plateau 
for small times. At such short times, the particles have 
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FIG. 3: (Color online) Correlation functions dj as defined in the text, and linear combinations thereof, comparing the 
experimental data (with error bars) with the numerical fit functions (without) for an optimized parameter set. The statistical 
error of the numerical data is smaller than the line width. Parts (a) - (f) have been obtained by modeling the observation 
volumes by Eq. [9l while for parts (g) and (h) we have assumed a Gaussian form (Eq. [G]). 
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essentially not moved at all. Hence the plateau indi- 
cates that a particle is able to send photons to both de- 
tectors from essentially the same position, or, in other 
words, that the effective observation volumes must over- 
lap quite substantially. This overlap effect then of course 
also shows up in the backward correlation function G21 
(see part (d)) at short times, with precisely the same 
plateau value. Therefore, such overlap effects essentially 
cancel out when considering the difference G12 — G21 in- 
stead (see part (f)), while of course they are strongly 
present in the mean (G12 -|- G2i)/2 (see part (e)). 

Obviously, the source of the overlap must be an effect 
of the optical imaging system, which is of course some- 
what complicated, due to the many components that are 
involved. However, beyond this general statement we 
have unfortunately so far been unable to trace down its 
precise physical origin, and therefore also been unable 
to construct a fully consistent model for the observation 
volumes. The simple models that we have considered in 
our present work are not fully adequate, meaning that 
they systematically underestimate the amount of over- 
lap, unless one assumes highly unphysical parameters, 
which would cause other aspects of the modeling to fail 
completely. It should be noted that similar overlap ef- 
fects are also present in standard double-beam FCCS 
[2^ : however, the underlying physics for that setup is 
slightly different, and the modeling used there cannot be 
simply transferred to our system. 

Fortunately, however, our best model for the observa- 
tion volumes is at least physical enough such that it can 
describe not only the autocorrelation functions (see part 
(b)) but also the overlap-corrected difference G12 — G21 
(part (f)) reasonably well, while still failing to describe 
the mean (G12 -|- G2i)/2 (part (e)). For this reason, our 
fitting procedure altogether takes into account the linear 
combinations (Gii-|-G22)/2 and G12 — G21, while deliber- 
ately discarding the data on (Gi2-l-G2i)/2 and Gn — G22- 
This is nicely borne out in Fig. |31 which shows not only 
the experimental data, but also the result of our theoret- 
ical modeling for optimized parameters. 

The fact that the success of the modeling depends cru- 
cially on an accurate description of the observation vol- 
umes is strongly underpinned by parts (g) and (h) of Fig. 
[31 The experimental data for (Gi2-l-G2i)/2 and G12— G21 
are again the same, but the theoretical model uses a dif- 
ferent functional form for the observation volumes, whose 
performance is obviously significantly poorer: Not only 
is the overlap plateau (part (g)) underestimated even 
more strongly than for the better model (part (c)), but 
also in the overlap-corrected function G12 — G21 (part 
(h)) are the deviations from the experimental data much 
more pronounced than for the better model (part (f)). It 
should also be noted that the autocorrelation functions 
are much less sensitive to these details; the autocorrela- 
tion curve for the poorer model (data not shown) fits the 
experiments as well as the better one (part (b)). 



IV. CORRELATION FUNCTIONS AND 
PARTICLE DYNAMICS 

A. Molecular Detection Efficiency 

The fluorescence particles pass consecutively through 
the two observation volumes Wi and W2 (Fig. [2]). The 
observation volume of each pinhole is given by the space- 
dependent molecular detection efficiency (MDE) func- 
tion. It depends on the excitation intensity profile Iz{z), 
and the collection efficiency of the objective plus detector 
system. In essence, the function Wi{r) denotes the prob- 
ability density for the event that a fluorescence photon 
emitted from a particle at position r will pass through 
pinhole 1 and reach detector 1. Similarly, W2{r) is the 
analogous function for pinhole 2. Since the intensity of 
the evanescent wave decays exponentially with a pene- 
tration depth dp (of order lOOnm), and the observation 
volumes are displaced with respect to one another by a 
distance Sx (roughly 800nm), we assume the functional 
form 

Wiir) = Wxy{x,y)d;^exp(^-£^ , (3) 
W2{r) = Wxy{x-sx,y)d;^cxp(^-£j , (4) 
where normalization of the probabihty densities imphes 

/oo poo 
dx / dyWxy {x, y) = 1. (5) 
-co J —OQ 

In general the function Wxy is given by the convolution 
of the pinhole image in the sample space with the point 
spread function (PSF) of the objective. However, one 
simple and widely used approximation, valid for pinholes 
equal or smaller than the Airy Unit of the system, as- 
sumes that Wxy is a Gaussian function p2, 32, iO,]: 

Wxyix, y) = ^ exp ("-2^^^^) ; (6) 

a typical value for the width that we obtain from fitting 
is Wo — 250nm. 

A substantially better description of Wxy can be ob- 
tained by considering the explicit form of the PSF 
[lllll^]. However, this form is described with complex 
mathematical equations and is often approximated by a 
squared Bessel function [13, EH, IH 

PSFxyO.{^-^)\ (7) 
where Ji denotes the first Bessel function and 



q^kNA ^x"^ + 2/2 = iVA ^a;2 + y"^ . (8) 
A 

Here A is the wavelength of the fiuorescent light (in 
our case GOOnm). The Bessel PSF implicitly assumes a 
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FIG. 4: (Color online) Comparison of the two normalized 
MDEs used in our study, for the optimized parameters of 
Fig. O using the natural unit system of the PCBPSF. 



paraxial approximation (i. e. small NA). While this 
assumption is probably not the best for confocal mi- 
croscopy, it is certainly more accurate than a simple 
Gaussian PSF ^E^. 

As mentioned above, in order to describe what a pin- 
hole sees one must calculate the convolution of the PSF of 
the objective with the pinhole image in the sample space. 
The geometrical image of the pinhole is simply obtained 
by dividing the physical size of the pinhole (physical ra- 
dius = 50iim) by the total magnification of the system 
(in our case « 333). This results in a radius RpH in the 
sample space of approximately ISOrim. Therefore the to- 
tal model MDE is given by ^ 



, , / kNA V f 



2„ fiJM 



where 



d'ro , 

\ro\<BpH \ 1 



(9) 



(10) 



The convolution integral is difficult to evaluate analyti- 
cally, but easy to calculate numerically. To this end, we 
use dimensionless length units in which the factor kN A 
is unity. In these dimensionless units, Rpu takes the 
value 2.3 for the parameters given above, which is the 
value we have used throughout our study. We call this 
function ^ the "pinhole-convoluted Bessel point spread 
function" (PCBPSF), which we calculated in dimension- 
less units once and for all, and stored as a table. During 
the actual data analysis, the transformation factor from 
dimensionless units to real units was used as a fit pa- 
rameter, in analogy to wq for the Gaussian model. It 
should be noted that the PCBPSF decays for large dis- 
tances like {x^ +2/^)"'^^^ J therefore providing much more 
overlap than the Gaussian model. 

In the present work, we have studied both models, 
the "Gaussian" model according to Eq. [HI as well as the 



PCBPSF model according to Eq. 121 The corresponding 
correlation curves have already been presented in Fig. |21 
The corresponding MDEs are shown in Fig. |4l One 
sees that the PCBPSF model puts much more statistical 
weight into the tail of the distribution than the Gaussian 
model. As already discussed above, we found the Gaus- 
sian model to perform less well than the PCBPSF model, 
since it underestimates the overlap even more severely 
than the latter. In what follows, we will present data 
always for the PCBPSF model, unless stated differently. 



B. Theory of Correlation Functions 

The dynamics of the tracer particles is described by 
the space- and time-dependent concentration (number of 
particles per unit volume) C{r,t), its fluctuation 

5C{r,t)^C{r,t)- {G) (11) 

and the concentration correlation function 

$(r,r',t)-(JC(r,t)JC(r',0)); (12) 

note that translational invariance applies only to time, 
but not to space, due to the presence of the flow and the 
surface. At time t = 0, this reduces to the static corre- 
lation function, for which we simply assume the function 
pertaining to an ideal gas: 



$(r,r',0) = {C)5{r-r'). 



(13) 



Note that this assumption implies that we consider the 
particles as point particles, with no interaction with the 
surface except impenetrability, and no interaction be- 
tween each other, due to dilution. 

As described in Ref. [l^l , the correlation functions are 
related to $ via 



/ / d^r(Pr'Wi{r')Wj{r)^{r,r' ,t) 



(C)2 {J (PrW^{r)) {J (PrWj[r)) 
= {Cy^ [ [ d^rd^r'W^{r')Wj{r)<^{r,r\t), 



(14) 



where in the second step we have taken into account the 
normalization of the Wi. Therefore, the obvious strategy 
for analyzing the experimental data is to (i) evaluate $ 
within a model for the particle dynamics, (ii) evaluate 
the integrals in Eq. [14] to obtain a theoretical prediction 
for Gij for a given set of parameters, (iii) compare the 
prediction with the data, and (iv) optimize the parame- 
ters. The normalizing prefactor (C)^^ is not known very 
accurately and will hence be treated as a fit parameter. 

The tracer particles undergo a diffusion process and 
move in an externally driven flow field v. Hence, 
we describe the concentration correlation function by a 
convection-diffusion equation of the form 

dt<i>ir,r',t) = i:>V^$(r,r',t)-V^-t;(r)$(r,r',t), (15) 
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which needs to be solved for z > 0, > with the initial 
condition Eq. [T2] and the no-flux boundary condition at 
the surface, 

9.$(r,r',i)L-=o=0, (16) 

which imposes that there is no diffusive current entering 
the solid. For reasons of simplicity, the hydrodynamic 
interactions with the surface are neglected, and hence the 
diffusive term is described only by an isotropic diffusion 
constant D. 

Since in the experiment the exponential decay length 
of the spatial detection volume normal to the surface is in 
the range of 100 — 200nm, while the channel size is three 
orders of magnitude larger, it is justified to assume the 
flow field to be approximately linear. For our geometry, 
this implies 

v{r)=^'i -{r + he,), (17) 

where Is is the slip length, 7 ~ dvx/dz is the constant 
shear rate, Bz denotes the unit vector in z-direction and 
^ = e-x ® e-z '^s the dimensionless rate-of-strain tensor. 

At this point, it is useful to re-define the coordinate 
system in such a way that the finite hard-sphere radius 
R of the tracer particles (roughly Inm) is taken into ac- 
count. We therefore identify z = no longer with the 
interface, but rather with the z coordinate of the par- 
ticle center at contact with the interface. In this new 
coordinate system, the flow field is given by 

v(r)=^'i ■{r + {h + R)ez), (18) 

i. e. we simply have to add the particle radius to the slip 
length. The functional form of the observation volumes 
Wi and W2 remains unchanged, since the z dependence is 
just an exponential decay, such that a shift in z direction 
just results in a constant prefaetor that can be absorbed 
in the overall normalization. Our method therefore does 
not yield a value for Is, but rather only for the combina- 
tion Is + R. 

As mentioned previously, for some special cases the 
convection-diffusion equation can be solved analytically, 
for example in the ease of uniform or linear flow in bulk, 
i. e. far away from surfaces [13, IH, El], or for pure 
diffusion close to the wall, but without any flow field 
[H,!!^. For our case, however, it is not easy, or even im- 
possible, to find such a solution. Therefore the aim of the 
next sections will be to construct a stochastic numerical 
method. Concerning the problems that were mentioned 
after Eq. [TH (i) and (ii) can be solved by Brownian Dy- 
namics, while problems (iii) and (iv) are tackled by a 
Monte Carlo algorithm in parameter space. 

V. SAMPLING ALGORITHM 

Brownian motion of particles under the influence of 
external driving is described by a Fokker-Planek equa- 
tion [47l - [5l| . which has exactly the same form as the 



convection-diffusion equation, Eq. [131 the only difference 
being that $ is replaced by the so-called "propagator" 
P(r,t|r',0), which is the conditional probability density 
for the particle motion r' — >■ r within the time t. P and 
$ describe the same physics and are actually identical 
except for a trivial normalization factor, $ = {C)P. We 
can therefore rewrite Eq. [14] as 

(C)G.,(t) (19) 




d\d^r'W^{r')Wj{r)P{r, t\r' , 0). 



As is well-known, the Fokker-Planek equation is equiv- 
alent to describing the particle dynamics in terms of a 
Langevin equation 

r{t) = v{r{t))+r^{t). (20) 

Here r{t) is the tracer velocity, v is the deterministic (ex- 
ternal) velocity imposed by the flow, while J7 is a stochas- 
tic Gaussian white noise term which describes the diffu- 
sion: 

Mt)) = 0, (21a) 
{i]o.{t')j^0{t)) = 2D5o.fi5{t' ~t). (21b) 

Here, a, /3 = x,y, z are Cartesian indices and da^ is the 
Kronecker delta. We solve this Langevin equation nu- 
merically by means of a simple Euler algorithm [49j with 
a finite time step At: 

r{t + At) = r{t) + Atv{r{t)) + V2DAtx, (22) 

where x = {Xx,XyiXz) is a vector of mutually indepen- 
dent random numbers with mean and variance 1. The 
boundary condition at the wall is taken into account by 
a simple reflection at 2: = 0, i. e. a particle that, after 
a certain time step, has entered the negative half-space 
2: < is subjected to 2; — >■ — z before the next propagation 
step is executed. 

Now, let us consider a computer experiment where, at 
time t = 0, we place a particle randomly in space, with 
probability density poir'), and then propagate it stochas- 
tically according to Eq. [211 The probability density for 
it reaching the position r after the time t is then given 
by 

Q{r,t)^ J £r'P{r,t\r',0)po{r'). (23) 

If we now consider an observable A, which is some func- 
tion of the particle's coordinate, A — A{r), and study the 
time evolution of its average, then this is obviously given 
by 

{A)it) = {Airm (24) 
= J (frA{r)Q{r,t) 




(fr(fr'A{r)P{r, t\r' , Q)po{r'). 
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Therefore, if we set po = Wi and A = Wj, then {A) 
is identical to the rescaled correlation function {C)Gij. 
In other words, we place the particle initially with prob- 
ability density Wi, then generate a stochastic trajectory 
via Eq. [521 and evaluate Wj for all times along that tra- 
jectory. This yields a function Wj{t) for that particular 
trajectory. This computer experiment is repeated often, 
and averaging Wj{t) over all trajectories yields directly 
a stochastic estimate for the (unnormalizcd) correlation 
function Gy . Of course, these estimates will have statisti- 
cal error bars, just as the experimental ones; however, we 
sample several hundred thousand trajectories, such that 
the numerical errors are substantially smaller than the 
experimental ones. In principle, the numerical data arc 
also subject to a systematic discretization error as a re- 
sult of the finite time step; however, by choosing a small 
value for A< we have made sure that this is still small 
compared to the statistical uncertainty. Note also that 
our approach implements an optimal importance sam- 
pling [521 with respect to the t = Q factor Wi, but not 
with respect to Wj. In practical terms, our straightfor- 
ward sampling scheme turned out to be absolutely ade- 
quate. 

The simulations were run using a "natural" unit sys- 
tem where length units are defined by setting dp to unity, 
while the time units are given by setting the diffusion 
constant D to unity. The time step was fixed in physi- 
cal units to a value of at most 2/is (it was dynamically 
adjusted in order to match the non-equidistant experi- 
mental observation times), which, for all parameters, is 
much smaller than unity in dimensionless units. Obvi- 
ously, this is small enough to represent the stochastic part 
of the Langevin update scheme with sufficient accuracy. 
For typical parameters {D = 35fj,m^/s, dp = 0.1/im, 
7 = 4 X 10^s~^), the dimensionless unit time corresponds 
to ~ 0.3ms, such that the resulting value for the dimen- 
sionless shear rate (~ 1.2) is of order unity as well. Since 
dp (or unity, in dimensionless units) defines the z range 



in which the statistically relevant part of the simulation 
takes place, we find that typical flow velocities in dimen- 
sionless units are also of order unity. This shows that the 
time step is also small enough for the deterministic part of 
the Langevin equation. We also see that the experiment 
is neither dominated by diffusion nor by convection, and 
therefore the analysis needs to take into account both. 

As a simple test case, we used our algorithm to calcu- 
late the autocorrelation function for vanishing flow and 
the Gaussian model for the observation volume, where 
an analytical solution is known jssl . |46| . In our dimen- 
sionless units, it is, up to a constant prefactor, given by 



G(^'(i) 



(25) 



(1 - 2t) exp (i) erfc 



Vt 




Figure [S] shows the analytic autocorrelation function with 
wq = 2 and its simulated counterpart, averaged over 10'^ 
independent trajectories, where a small time step of At = 
lO"'^ (in dimensionless units) was used. In Fig. [5] the 
deviation of the simulated data {G^^'' ) from the analytic 
expression is shown, 



error (t) ^ G^'Ht) - G'^^'^t). 



(26) 



Clearly, the numerical solution converges to the analyt- 
ical result when the number of trajectories is increased, 
as it should be. 



VI. SUBTRACTION SCHEME 

At this point, it is worthwhile to reconsider the sub- 
traction procedure introduced in Sec. IIIIl To this end, we 
assume that the true functions W; differ somewhat from 



the model functions, which we will denote by Wj 
is most easily parameterized by the ansatz 



(m) 



Wi = (1 -e)W;^'"' +eW, 



This 



(27) 
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where Wi,!^^"^ and Wi are all normalized to unity, while 
e is a (hopefully) small parameter. For the purposes of 
the present analysis, we also assume that the Brownian 
Dynamics model is a faithful and correct description of 
the true dynamics, i. e. that the difference between Wi 

and M^j^™'' is the only reason for a systematic deviation 
between simulation and experiment. 

Inserting Eq. [57] into Eq. [121 we thus find 

{C)G,j{t) (28) 

= {l-ef I f d\d\'wt'\r')wl'''\r)P{r,t\r',0) 



+ e(l-£) / I d^rd^r'W^^"'\r')Wj{r)P{r,t\r',0) 
+ £(l-£) J y dVdVVF,(r')W^j™'(r)F(r,i|r',0) 
+ [ [ d^rd^r'W^{r')Wj{r)P{r,t\r\0). 



Since we treat (C) as an adjustable parameter, it makes 
sense to view the first term (including the prefactor 
(1 — e)^) as the theoretical model for the correlation func- 
tion, (C) G[^\t). For the deviation between experiment 
and theory we then obtain, neglecting all terms of order 

£^ 



K, 



'\C) ( 



(29) 



d^rd^r'wl""^ {r')Wj (r)P(r, t\r' , 0) 



+ j j d'^rd^r'W,{r')wl"'\r)P{r,t\r' ,0), 
and for its antisymmetric part 



(30) 



d\d\'[W^^"'\r')W,{r) 



W,*"^(r)VK,-(r')]P(r,t|r',0) 
d^rd^r'[wj"'\r')mir) 

-W^"'\r)m{r')]P{r,t\r',0). 

The terms in square brackets are antisymmetric under 
the exchange r r' , and hence P can be replaced by its 
antisymmetric part 

Pa(r, r', t) = P(r, t\r\ 0) - P{r' , t\r, 0). (31) 

Exchanging the arguments in the second terms within 
the square brackets then yields 



d^rd^r'W^'''\r')Wj{r)Pa{r,r',t) 
dVdViyf"^(r')W'»(r)Pa(r,r',t). 



(32) 



This is clearly a nonzero contribution. In other words, 
the subtraction scheme (i. e. studying G12 — G21 instead 
of G12) does not provide a consistent cancellation pro- 
cedure such that the first-order deviation would vanish. 
However, in practical terms the deviation is much smaller 
than for the original data (G12 and G21), for which Eq.[25| 
applies. To some extent, this is so because the error 
is the difference of two terms, but mostly it is due to 
the fact that not the full propagator P contributes, but 
rather only its antisymmetric part Pa- For short times 
the dynamics is dominated by diffusion, i. c. P is es- 
sentially symmetric, or Pa ~ 0. At late times, we again 
expect Pa to become quite small (exponentially damped, 
see Eq. IM)) . although we have no rigorous proof for this. 
Therefore one should expect that the strongest deviation 
occurs at intermediate times where Pa is maximum. This 
time scale is not given by the optical geometry but rather 
by the dynamics; dimensional analysis then tells us that 
this time must be of order D/v^ . For typical parameters 
of our experiment [D = 35^m^/s, w = 4 x IO^/xto/s) 
we obtain a value of roughly 0.2ms, which fits quite well 
to the observations one can make in Fig. [3l part (h). 
At such times, we expect that the main contribution to 
A'i2 — K21 comes form the first term of Eq. [321 (down- 
stream vs. upstream correlation) and that Pa is positive 
for most of the relevant arguments. Therefore, one should 
expect that the experimental data should lie systemat- 
ically above the theoretical predictions, which is indeed 
the case. Our expectations concerning the behavior of Pa 
come from studying the simple case of one-dimensional 
diffusion with constant drift without boundary condi- 
tions; here one has 



P{x,t\x',()) 



1 



cxp 



[x — x' — vty 

4Dt 



and 



Pa {x^ X , ^) 



exp 



{x-x'f 



exp 



vH 



sinh 



4:Dt 

{x — x')v 
2D 



(33) 



(34) 



VII. STATISTICAL DATA ANALYSIS 

A. Monte Carlo Algorithm 

For the model that we consider in the present pa- 
per, the space of fit parameters is (in principle) seven- 
dimensional. We have three lengths that define the ge- 
ometry of the optical setup, dp, Sx, and wq (Gaussian 
model) or {kNA)~^ (diffraction model). Three further 
parameters define the properties of the fiow and the diffu- 
sive dynamics of the tracers; these are the diffusion con- 
stant D, the shear rate 7, and the slip length plus particle 
radius Is + R. Finally, there is the concentration of tracer 
particles (G) , which serves as a global normalization con- 
stant. The functions to be fitted are (Gn + G22)/2 and 
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G12—G21 ■ However, we have seen in Scc. lVII that the non- 
idealities in modeling the observation volumes do have 
an effect on the normalizations, and therefore we allowed 
one separate normalization constant (C) for each of the 
curves ((C) a for the autocorrelation and (C)c for the 
cross-correlation), in order to partly compensate for these 
non-idealities. Therefore, our parameter space is finally 
eight-dimensional. The strategy that we develop in the 
present section aims at adjusting all parameters simulta- 
neously in order to obtain optimum fits. For the further 
development, it will be useful to combine all the parame- 
ters into one vector n. Furthermore, for each parameter 
we can, from various physical considerations, define an 
interval within which it is allowed to vary (because val- 
ues outside that interval would be highly unreasonable 
or outright unphysical). This means that we restrict the 
consideration to a finite eight-dimensional box 51n in pa- 
rameter space. 

A central ingredient of our approach is the fact that 
both the experimental data and the simulation results 
have been obtained with good statistical accuracy (~ 
2.5 X 10^ trajectories for the simulations, 40 independent 
measurements for the experiments). This does not only 
allow us to obtain rather small statistical error bars, but 
also (even more importantly) to rely on the asymptotics 
of the Central Limit Theorem, i. e. to assume Gaussian 
statistics throughout. For both correlation curves and 
each of the considered times, we have both an experimen- 
tal data point Ei and a simulated data point Si, where 
the index i simply enumerates the data points. Both Ei 
and Si can be considered as Gaussian random variables 
with variances a'^ ^ and (t| ^ , respectively. Then 

A, = (35) 

is again a Gaussian random variable, whose variance is 
simply unity, 

(a^) - (A,)' = 1. (36) 

Therefore. A^ is, in principle, a perfect variable to mea- 
sure the deviation between simulation and experiment. 
Unfortunately, however, the parameters <Ts,i and asA are 
not known. What is rather known are their estimators 
ss^i and SE,i, as they are obtained from standard analysis 
to calculate error bars. Therefore, we rather consider 

A, = ^'^^^ . (37) 

The statistical properties of this variable, however, are 
in the general case unknown (53l |. It is only in the case 
of rather good statistics (as we have realized it) that we 
can ignore the difference between a and s, and simply 
assume that A^ is indeed a Gaussian variable with unit 
variance. It is at this point where the statistical quality 
of the data clearly becomes important. 



If M is the total number of data points, then 

M 

n^-Y^A^ (38) 

1=1 

is obviously a quantity that measures rather well the de- 
viation between experiment and simulation. In princi- 
ple, the task is to pick the parameter vector 11 in such 
a way that T-l is minimized. We have deliberately cho- 
sen the symbol T-l in order to point out the analogy to 
the problem of finding the ground state of a statistical- 
mechanical Hamiltonian. In case of a perfect fit, we have 
{S^) = {E,) or (A,) = 0, implying {H) = M/2. In the 
standard nomenclature of fitting problems, 2H is called 
"chi squared". We also introduce ^ = 2H/M, which we 
will call the "goodness of simulation" (standard nomen- 
clature: "chi squared per degree of freedom" ) . 

For optimizing 11, we obviously need to consider Ji as 
a function of 11. In this context, it turns out that it is 
important to be able to consider it as a function of only 
n, and to make sure that this dependence is smooth. 
For this reason, we use the same number of trajectories 
when going from one parameter set to another one, and 
use exactly the same set of random numbers to generate 
the trajectories. In other words, the trajectories differ 
only due to the fact that the parameters were changed. 
Therefore, both Si and ss,i are smooth functions of the 
parameters, and H is as well. 

In order to find the optimum parameter set, one could, 
in principle, construct a regular grid in fljj and then 
evaluate Tl for every grid point. However, for high- 
dimensional spaces (and eight should in this context be 
viewed as already a fairly large number, in particular 
when taking into account that it is bound to increase 
further as soon as more refined models are studied), it is 
usually more efficient to scan the space by an importance- 
sampling Monte Carlo procedure based upon a Markov 
chain [52|. Applying the standard Metropolis scheme 
[5^ . we thus arrive at the following algorithm: 

1. Choose some start vector 11. This should be a rea- 
sonable set of parameters, perhaps pre-optimizcd 
by simple visual fitting. 

2. From the previous set of parameters, generate a 
trial set via II' = 11 + AH, where AH is a random 
vector chosen from a uniform distribution from a 
small sub-box aligned with fin- 

3. If the new vector is not within 57n, reject the trial 
set and go to step 2. 

4. Otherwise, calculate both Peg(n') and Peq{Il), as 
well as the Metropolis function 

m = min (l, P.^iji) j P,q{W)) , (39) 

where Peg is the "equilibrium" probability density 
of n, i. e. the desired probability density towards 
which the Markov chain converges (more about this 
below) . 
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5. Accept the trial move with probability m (reject it 
with probability 1 — m) , count either the accepted 
or the old set as a new set in the Markov chain, and 
go to step 2. 

6. After relaxation into equilibrium, sample desired 
properties of the distribution of 11, like mean val- 
ues, variances, covariances, etc., by simple arith- 
metic means over the parameter sets that have been 
generated by the Markov chain. This allows the es- 
timation of not only the physical parameters, but 
at the same time also of their statistical error bars. 

The scheme is defined as soon as Peg is specified. Now, 
from the considerations above, we know that in case of 
a perfect fit the variables are independent Gaussians 
with zero mean and unit variance. This implies (ignor- 
ing constant prefactors which anyway cancel out in the 
Metropolis function) 

Peg cx Y[exp(-^AA 

= exp(~H), (40) 

which makes the interpretation in terms of statistical me- 
chanics obvious. Clearly, this form for P^q is the only 
reasonable choice for implementing the Monte Carlo al- 
gorithm. After relaxation into equilibrium, one should 
observe a ^ value of roughly unity, while larger numbers 
indicate a non-perfect fit (even after exhaustive Monte 
Carlo search), and thus deficiencies in the theoretical 
model. One should also be aware that the equilibrium 
fluctuations of ^ are expected to be quite small, since 
^ is the arithmetic mean of a fairly large number (M, 
the number of experimental data points) of independent 
random variables. 

In practice, we adjusted AH in order to obtain a fairly 
large acceptance rate of roughly 0.6. . .0.8. The Monte 
Carlo algorithm was then run for more than 3 x 10^ steps, 
each step involving the generation of roughly 2.5 x 10^ 
trajectories. The simulation was run on 512 nodes (2048 
processes) of the IBM Blue Gene-P at Rechenzentrum 
Garching, where each process generated 123 trajectories. 
On this machine, one Monte Carlo run took roughly one 
day to complete. It turned out that discarding the first 
5 X 10^ configurations was sufficient to obtain data in 
equilibrium conditions, where the mean values of the pa- 
rameters and their standard deviations were calculated. 
It should be noted that the equilibrium fluctuations of 
the parameters tell us the typical range in which they 
can still be viewed as compatible with the experiments. 
Therefore these fluctuations are the appropriate measure 
to quantify the experimental error bars, while calculating 
a standard error of mean (or a similar quantity) would 
not be appropriate and severely underestimate the errors. 
Finally, it should be noted that the approach allows in 



principle to analyze the mutual dependence of the param- 
eters as well, by sampling the corresponding covariances; 
this was however not done in the present study. 

B. Scale Invariance 

As noted before, the correlation functions depend on 
the average concentration (C), the diffusion constant D, 
the shear rate 7, and various lengths, which we denote 
by Simple dimensional analysis shows that for any 
scale factor a the scaling relation 

G,,{t,{C),D,^,{k}) (41) 
= G,,{t,a^C),D/a\^,{l,/a}) 

holds. The "Hamiltonian" of the previous subsection is of 
course subject to the same scale invariance. This means 
that for each point in parameter space there is a whole 
"iso-line" in parameter space that fits the data just as 
well as the original point. Therefore, in order to im- 
prove the MC sampling, we generated such an iso-line 
for each point in parameter space that was produced by 
the Markov chain of the previous subsection. Of course, 
the iso-lines were confined to the region of the overall 
parameter box. It turned out that our Markov chains 
were still so short that this improvement was not com- 
pletely superfluous (as it would be in the limit of very 
long chains). In other words, taking the invariance into 
account helped us to avoid underestimating the errors. 

In practice, this was done as follows: Assuming that 
the most accurate input parameters are the penetra- 
tion depth dp = 100 ± lOnm, the diffusion constant 
-D = 36 ± b^m?/s and the separation distance — 
800 ± SOnm, we calculate for every data point a mini- 
mum and a maximum scaling factor a, such that we ob- 

, . j(min) _i 7 i(rnax) (min) _i (max) 

tam dp < a dp < dp , si < a < Sx 
and d(™") < a-^D < for all a in (a„„„, a„a:r)- 

This provides us with additional data points in parame- 
ter space that are added to the statistics. 

C. Sample-to-Sample Fluctuations 

It should be noted that the parameters found by the 
procedure outlined above are optimized for a specific set 
of random numbers used to generate the trajectories. 
Therefore, one must expect that one obtains different 
results when changing the set of random numbers. In 
our statistical-mechanical picture, we may view the set 
of random numbers as random "cou plin g constants" of 
a disordered system like a spin glass [5j] , where the dis- 
order is weak since the number of trajectories is large. 
For disordered systems, the phenomenon of "sample- 
to-sample" fluctuations is well-known, and it should be 
taken into account. We have therefore run one test where 
we applied the same analysis to three different random 
number sequences. Indeed we found (see Tab. |T| that 
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seed 


42 


4711 


2409 




av 


a 


av 


a 


av 


a 




17.83 


1.92 


17.81 


1.96 


17.94 


1.94 


{C)c[t^m-'''] 


16.53 


1.80 


16.45 


1.82 


16.83 


1.83 


dp [nm] 


95.81 


3.50 


96.02 


3.58 


95.78 


3.52 


[kNA)-^[nm] 


68.70 


2.60 


68.58 


2.61 


69.81 


2.69 


Sx [nm] 


781.94 


27.70 


779.54 


28.28 


793.22 


28.24 




36.59 


2.56 


36.47 


2.62 


36.63 


2.57 


Is + -R[nm] 


12.80 


1.10 


11.62 


0.90 


15.14 


1.21 




1.441 


0.018 


1.277 


0.016 


1.718 


0.017 


acceptance rate 


82.0% 




82.0% 




82.2% 




No. of MC steps 


609410 




609590 




610510 





TABLE I: Averaged values (av) and standard deviations (a) 
calculated from MC simulations with fixed 7 — 3800s~^, but 
different start values ( "seeds" ) for the random number gener- 
ator. 
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FIG. 7: Goodness of simulation ^ as function of the number 
of Monte Carlo steps for 7 = 3800s^^. 



sample-to-sample fluctuations are observable, and some- 
what larger than the errors obtained from simple MC, 
while still being of the same order of magnitude. A con- 
servative error estimate should therefore take these fluc- 
tuations into account, by multiplying the error estimates 
from plain MC by, say, a factor of three. In what fol- 
lows, we will only report the simple MC estimates for 
the errors. 



VIII. RESULTS 



The experiments were performed with a penetration 
depth of the evanescent wave of dp ~ lOOnm, the lat- 
eral size of the observation volumes (within the Gaussian 
model) was wq ~ 250nm and their center-to-center sepa- 
ration was ~ 800nm. Furthermore, the diffusion con- 
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FIG. 8: Slip length plus particle radius as function of the 
number of Monte Carlo steps for 7 — 3800s^^. 



slant of the tracers is known to be roughly D ~ 36^m^/s 
as measured by dynamic light scattering. The shear rate 
was determined from an independent measurement using 
single-focus confocal FCS [2^ where the entire flow 
profile across the microchannel was mapped out. Alter- 
natively, one might also use double-focus confocal FCCS 
[13, Hyj- From this measurement, we obtained a shear 
rate at the bottom channel wall of 7 = 3854 ± 32s~^. 
More details on this issue and some theoretical back- 
ground are presented in the appendix. Nevertheless, 
we took a conservative approach and allowed the shear 
rate to vary between 3500s~^ and 4000s~^. Finally, 
we expected the slip length to be not more than a few 
nanometers, but we nevertheless allowed it to vary up to 
~ lOOnm. These estimates allowed us to start the Monte 
Carlo procedure with good input values. 

We then observed the Monte Carlo simulation to sys- 
tematically drift to smaller and smaller values of 7, until 
finally "getting stuck" at the imposed lower boundary, 
7 = 3500s~^. What we mean by this term is a behavior 
where fiuctuations near 3500s~^ still occur, but in such 
a way that 3500s~^ is the most probable value, while 
smaller values only do not occur because we do not allow 
them. Since we know experimentally that 7 = 3500s~^ is 
clearly unacceptable, this behavior again indicates that 
the theoretical model is not completely sufficient to de- 
scribe the experimental data (sec also the discussion in 
Sees, [m] and Ell)- 

We therefore decided to keep 7 fixed during a Monte 
Carlo run, and rather vary it systematically in the given 
range. For none of the parameters were we able to obtain 
a better goodness of simulation than ^ ~ 1.25, which is 
still a bit too large, i. e. indicates a non-perfect fit (al- 
though the data on (Gi2+G2i)/2 have been discarded al- 
ready) . The convergence behavior of the method is shown 
in Fig. [7J where we plot ^ as a function of the number 
of Monte Carlo iterations. For the Gaussian model, the 
best ^ value that we could obtain was ^ ~ 2.5, which is 
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FIG. 9: Averaged slip length as function of the shear rate, 
calculated from the Monte Carlo results. 



substantially worse. 

With these caveats in mind, we may proceed to study 
the parameter values that the Monte Carlo procedure 
yields. Obviously, the most interesting one is the slip 
length Is , or the sum ls + R (recall that the method docs 
not provide an independent estimate for these parame- 
ters, but only for their sum). Figure [8] presents data on 
the evolution of Zs + i? during the Monte Carlo process 
for 7 = 3800s~^; Is + R is thus seen to fluctuate between 
roughly lOnm and Ibnm, which is, within the limitations 
of the model, the statistical experimental uncertainty of 
this quantity. The mean and standard deviation of Is + R 
is shown in Fig. 12] as a function of shear rate, which are 
thus clearly seen to not be independent. Since we know 
7 much more accurately than the range plotted in Fig. |9l 
we see that in principle a fairly accurate determination 
of Is is possible, if the underlying theoretical model is de- 
tailed enough to fully describe the physics. One should 
note that the particle size R (more precisely, its hydrody- 
namic radius) is roughly 7nm; taking this into accunt as 
well, we find a value that is clearly smaller than lOnm. 
One should also note that for the Gaussian model we 
obtained a very similar curve; however, here the Ig + R 
values are systematically smaller by roughly 5nm. This 
again highlights the importance of having an accurate 
model for the MDE. 

The other results obtained from our MC fits are re- 
ported in Tab. [TTl 

Clearly, the Is values of Fig. [2] could only be viewed 
as definitive experimental results on Is if the agreement 
between experiment and model were perfect, with ^ ~ 1, 
and a good fit of all correlation functions. The reasons 
for the observed deviations are not completely clear; how- 
ever, all our findings hint very strongly to deficiencies in 
the description of the observation volumes, i. e. too inac- 
curate modeling of the detailed optical phenomena that 
finally give rise to the shape of these functions. Neverthe- 



7[s 



dp [nm] 

{kNA)-^[nm] 
Sx [nm] 
D[fim^/s] 
la + R[nm] 

acceptance rate 
No. of MC steps 



7[s 



{c)A[^^m■-^] 

{C)clt^m-'] 
dp [nm] 

{kNA)-^[nm] 
Sx [nm] 

Is + R[nm] 

acceptance rate 
No. of MC steps 



3500 



17.82 1.90 
19.94 1.82 

95.83 3.46 

68.84 2.53 
774.39 26.59 
36.71 2.49 
21.92 1.31 
1.377 0.016 
82.1% 
608090 



3800 



17.83 1.92 

16.53 1.80 

95.81 3.50 

68.70 2.60 

781.94 27.70 

36.59 2.56 

12.80 1.10 

1.441 0.018 

82.0% 

609410 



3600 



17.91 1.88 

17.00 1.80 

95.66 3.42 

69.10 2.55 

777.98 26.63 

36.76 2.48 

19.16 1.24 

1.40 0.016 

82.1% 

611290 



3900 



17.67 2.00 

16.21 1.86 

96.14 3.70 

68.42 2.71 

783.75 28.76 

36.46 2.64 

9.92 1.12 

1.464 0.019 

81.9% 

610330 



3700 



17.93 1.89 

16.82 1.78 

95.64 3.44 

69.01 2.56 

780.53 26.97 

36.72 2.51 

15.98 1.12 

1.420 0.017 

82.1% 

610160 



4000 



17.96 1.89 

16.41 1.74 

95.58 3.41 

69.01 2.56 

789.77 27.35 

36.71 2.51 

7.88 0.99 

1.477 0.016 

82.0% 

610440 



TABLE II: Averaged values (av) and standard deviations (a) 
calculated from MC simulations with various shear rates. 



less, one should also bear in mind that the dynamic model 
is also rather simple, neglecting both hydrodynamic and 
residual electrostatic interactions with the wall. While 
one must expect that further refinements of the model 
will change both the Is values as well as their error bars, 
we believe that it is not probable that such a change 
would be huge. Given all the various systematic uncer- 
tainties of the modeling, we would, in view of our data, 
not exclude a vanishing slip length, while we consider a 
value substantially larger than, say, 15nm as fairly un- 
likely. 

Let us conclude this section by a few more remarks 
concerning our choice of parameters and the systematic 
errors of the method. From the setup it is clear that there 
are three parameters that can be varied experimentally 
fairly easily — these are the shear rate 7, the penetra- 
tion depth dp, and the effective pinhole-pinhole distance 
in sample space, Sx- The choice of parameters was gov- 
erned by various experimental considerations, which we 
will attempt to explain in what follows. 

It is clear that one wants a fairly large shear rate 7, 
in order to ensure that the signal has a sizeable contri- 
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bution from flow effects. In practice, liowever, increasing 
7 further by a substantial amount is limited by experi- 
mental constraints, such as channel construction, beaker 
elevation, etc. Furthermore, the choice of dp is subject 
to similarly severe experimental constraints: Increasing 
dp substantially would mean that we would approach the 
limit angle of total reflection closely, which would result 
in a very inaccurate a priori estimate of dp. On the other 
hand, an even smaller penetration depth value would be 
too close to the limits of the capabilities of the objec- 
tive, resulting in possible optical distortion effects which 
we would like to avoid. Finally, the choice of was gov- 
erned by our early attempts to suppress overlap effects by 
simply picking a fairly large value, such that the overlap 
integral is small. There are however two problems about 
such an idea. Firstly, a large value of decreases not 
only the overlap, but also the cross-correlation function 
as a whole [34l | , such that it ultimately becomes impossi- 
ble to sample the data with sufficient statistical accuracy 
on the time scale on which we can confidently keep the 
experimental conditions stable. Our value of should 
therefore be viewed as limited by such considerations. 
However, secondly, and more importantly, we realized 
in the course of our analysis that the overlap issue is 
not a problem of an unintelligent choice of parameters at 
all, but rather of our insufficient theoretical modeling of 
the MDE functions. As we have seen above, our results 
for the slip length depend rather sensitively and fairly 
substantially on the choice of the MDE function (up to 
nearly a factor of two). In our opinion, there is no rea- 
son to assume that this dependence would go away if we 
had picked parameters in a regime of small or vanishing 
overlap. From this perspective, we view the overlap es- 
sentially as a blessing, since it shows us where the main 
source of systematic error is probably located. 



In the light of these remarks, it would of course be 
interesting to systematically investigate the influence of 
the parameters dp and Sx on our results. In terms of 
the correlation functions as such, this has been done in 
Rcf. 



34| . and we refer the interested reader to that pa- 



per. However, doing the full analysis for a whole host 
of parameters would imply a very substantial amount of 
work, since all the experimental curves would have to be 
re-sampled again, in order to meet the rather stringent 
requirement of statistical accuracy that is built into our 
approach. We have hence not attempted to do this, but 
rather believe that it will be more fruitful to concentrate 
the efforts of future work on attempts to improve the the- 
oretical MDE modelling, even if that will be challenging. 
As far as the slip length is concerned, one must of course 
expect that the fitted value will depend on parameters 
such as dp and Sx, but only to the extent that this re- 
flects the systematic error — if the physics were modeled 
perfectly correctly, we would of course always obtain the 
same value. 



IX. CONCLUSIONS 

The results from the previous sections demonstrate 
that the method of TIR-FCCS in combination with the 
presented Brownian Dynamics and Monte Carlo based 
data analysis is in principle a very powerful tool for the 
analysis of hydrodynamic effects near solid-liquid inter- 
faces. Already within the investigated simple model of 
the present paper, we can conclude that the slip length 
at our hydrophilic surface is not more than lOnm. It was 
only the data processing via the Brownian Dynamics / 
Monte Carlo analysis that was able to demonstrate how 
highly sensitive and accurate TIR-FCCS is. 

The computational method has the advantage to be 
easily extensible to include more complex effects. For 
example, the hydrodynamic interactions of the particles 
with the wall would cause an anisotropy in the diffu- 
sion tensor [HH and a z dependence, electrostatic inter- 
actions would give rise to an additional force term in the 
Langevin equation, while polydispersity could be investi- 
gated by randomizing the particle size and the diffusion 
properties according to a given distribution. While these 
contributions are expected to yield a further improve- 
ment of the method, this was not attempted here, and 
is rather left to future investigations. However, we have 
also identified the inaccuracies in modeling the observa- 
tion volumes as (most probably) the main bottleneck in 
finding good agreement between theory and experiment, 
i. e., at present, as the main source of systematic errors, 
which makes it difficult to find a fully reliable error bound 
on the value of the slip length. 

Conversely, the problem of dealing with statistical er- 
rors can be considered as solved. For an extensive data 
analysis, as it has been presented here, one may need 
a supercomputer in order to obtain highly accurate re- 
sults in fairly short time. Nevertheless, the method will 
yield meaningful results even if confined to just a single 
modern desktop computer. Given the moderate amount 
of computer time on a high-performance machine, one 
should expect that quite accurate data should be ob- 
tainable within reasonable times by making use of the 
powerful newly emerging GPGPU cards. 

In our opinion, the presented method is a conceptually 
simple and widely applicable approach to process TIR- 
FCCS data, that is clearly only limited by inaccurate 
modeling. We believe that it has the potential to be- 
come the standard and general tool to process such data, 
in particular as soon as the optics is understood in better 
detail. The principle as such is applicable to all kinds of 
correlation techniques, such as FCS/TIR-FCS etc., and 
we think it is the method of choice whenever one inves- 
tigates a system whose complexity is beyond analytical 
treatment. 
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Appendix A: Solution of the Stokes Equation in a 
Rectangular Channel 

The flow profile throughout the height of the mi- 
crochanncl was measured by single-focus FCS under the 
same conditions as the TIR-FCCS experiments; the re- 
sult is shown in Fig. (TU] From a fit via a Poiseuille profile 
(solid line), we obtained an independent estimator for the 
shear rate near the wall, 7 = 3854 ± 32s~^. 

The purpose of this appendix is to analyze the theo- 
retical background of this fit in some more detail. For 
a pure Poiseuille flow, i. e. a simple parabolic profile, it 
is clear that the shear rate at the surface docs not de- 
pend on the slip length /j, because in this case a finite 
Is value simply shifts the profile by a constant amount. 
Therefore, in this case Is is indeed irrelevant for the fit. A 
short discussion on such issues is also found in Ref. po} . 
and experimentally [13, it is also known that typi- 
cally the shift is so small that a finite slip length is hard 
to detect by direct measurements of the profile. How- 
ever, from a theoretical and quantitative point of view 
it is not quite clear how well it is justified to assume a 
strictly parabolic profile, i. e. to assume that the flow 
extends infinitely in y direction — in our experiments, 
Ly/Lz ~ 40, which is large but not infinite. For finite 
values of Ly/Lz, the profile is somewhat distorted, and 
if this distortion is sufficiently large, then also a possible 
effect of Is should be taken into account. These questions 
can be easily answered by solving the flow problem in a 
rectangular channel in the presence of slip exactly, and 
this shall be done in what follows. The result of this anal- 
ysis will be that for our conditions the distortion of the 
profile is indeed completely negligible, and that therefore 



Is needs not be taken into account either. 
We start by considering the Stokes equation 



7? 



92 32 



(Al) 



in a rectangular channel with dimensions [—Ly/2, Ly/2] x 
[—Lz/2,Lz/2] in the yz-plane, as in the experiment. 
Here, 77 is the viscosity of the liquid and / is the driving 
force density or pressure gradient acting on the liquid in 
x-direction. We assume that all surfaces have the same 
slip length. 

For the case of a no-slip boundary condition, the solu- 
tion has been given in the textbook of Spurk and Akscl 
[56| . however in a form that docs not explicitly spell out 
the symmetry under exchange of y and z. Here we give 
the solution in a form that shows that symmetry, and 
generalize it to the case of a nonvanishing slip length Ig. 

Using the methods and notation of quantum mechan- 
ics, and allowing for some minor amount of numerics to 
evaluate a series, the solution is simple and straightfor- 
ward. We identify a function /(j/, z) with a vector |/) in 
a Hilbert space, and define the scalar product as 



U\9) 



+Ly/2 j- + L^/2 

dy / dzf*{y,z)g{y,z). 

-Ly/2 J-LJ2 



Defining a "Hamilton operator" via 
the Stokes equation is written as 



n\v^) 



ID 



Obviously, the functions 

\ky,kz) = N{ky, kz) cos{kyy) cos{kzz) 



(A2) 



(A3) 



(A4) 



(A5) 



with ky > 0, kz > and 



Niky,kz) 



Ly ^ SmjkyLy) 

Lz_ svaikzLz) 
2 2K 



-1-1/2 



'1/2 



are normalized eigenfunctions of Ti. , 



n\ky,kz)=j{kl + kl)\ky,kz) 



with 



{ky,kz\qy,qz) = 4^ 



,<ly 



(A6) 



(A7) 



(A8) 



The eigenfunctions must satisfy the boundary condi- 
tions, and hence the discrete wave numbers ky, kz must 
be the solutions of the transcendental equations 



/ k 

is^y 



^skz 



cot 
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2 
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17 



0.125- 



0.1 

"•"■I s-0.075- 



0.05 



0.025- 




-0.5 -0.4 -0.3 -0.2 -0.1 0.1 0.2 0.3 0.4 0.5 



FIG. 11: (Color online) One-dimensional cut of the flow pro- 
file at y = for no-slip boundary conditions and several values 
of Ly. 




allows us to rc- write Eq. IA6l as 



N{ky,k,) 



^ + Z,sin2 ( ky^ 



— -l-issm — 



-1/2 



-1/2 



(AlO) 



Since the set of eigenfunctions is complete, the spectral 
representations of % and "H^^ are given by 

^ = 7 E (^y + ^^') 1^^' ^^y^ ^-1 ' (All) 



\ky^ kz) {ky, kz I , (A12) 



resulting in the solution 



(A13) 



{ky, kz |1) \ky, kz) 
4 



N{ky,k,y 



ky kz 



X sin ( ky-Yj ^"^^ 



cos(fcj^y) cos(fc2z). 



Figure [TT] shows the resulting flow profile at j/ = (in 
the center of the channel) as a function of z, for vanishing 
slip length and various width-to- height ratios Ly/Lz of 
the channel. One sees that the convergence to the asymp- 
totic Poiseuille profile wp is indeed extremely rapid. The 
deviation, defined via 



FIG. 12: Averaged deviation between a one-dimensional cut 
of the fiow profile at y = for no-slip boundary conditions 
and the Poiseuille solution, as function of the width-height 
ratio of the channel. 



which, in the general case, can be found numerically. In 
the no-slip case, the solutions are simply given by ky = 
IT / Ly,?>T: / Ly, . . . and analogously for kz- Equation IA9I 



relative error : 



Lz 



dz 



vp{z) 



(A14) 



is displayed as a function of Ly/Lz in Fig. [T^ The 
rate of convergence is apparently exponential, and for the 
experimental value Ly/Lz = 40 the deviation is seen to 
be much smaller than the resolution of the measurements. 
Therefore, the assumption of a parabolic profile is indeed 
justified. 
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